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I. INTRODUCTION 

Simulations in lattice QCD have advanced remarkably in the past couple of years reaching within 100 MeV of 
the physical pion mass. This progress is due to theoretical improvements in defining the theory on the lattice and 
to algorithmic improvements that give a better scaling behavior as the quark mass decreases. These developments, 
combined with the tremendous increase in computational power, have made ah initio calculations of key observables on 
hadron structure in the chiral regime feasible enabling comparison with experiment. The hadron mass spectrum [1, 2] 
illustrates the good quality of lattice results that can be obtained. The agreement with experiment is a validation of 
the lattice approach, and justifies the computation of hadron observables beyond hadron masses, such as form factors 
and parton distribution functions. Both form factors and parton distribution functions can be obtained from the 
so-called generalized parton distributions (GPDs) in certain limiting cases. GPDs provide detailed information on the 
internal structure of hadron in terms of both the longitudinal momentum fraction and the total momentum transfer 
squared. Beyond the information that the form factors yield, such as size, magnetization and shape, GPDs encode 
additional information, relevant for experimental investigations, such as the decomposition of the total hadron spin 
into angular momentum and spin carried by quarks and gluons. 

GPDs are single particle matrix elements of the light-cone operator [3, 4], 

where n is a light-cone vector, and V denotes a path-ordering of the gauge fields in the exponential. Such matrix 
elements cannot be calculated directly in lattice QCD. However, 0{x) can be expanded in terms of local twist-two 
operators 

Q{,{i^^ti2-tin} ^ ^/pfMij D^'^■■■iD (2) 

where D = -^{D — D ) and {/ii, • • • , denotes symmetrization of indices and subtraction of traces. In this work 

we focus on the Dirac structures F = 7^^, 75 7'' and 7^ cr^"^ (cr'"^ = [7^, 7'^]/2), which are referred to as vector Oy{x), 

axial- vector 0^{x) and tensor 0^{x) operators, respectively. In lattice QCD we consider matrix elements of such 
bilinear operators. A number of lattice groups are producing results on nucleon form factors and first moments of 
structure functions closer to the physical regime both in terms of pion mass as well as in terms of the continuum 
limit [5-11]. While experiments are able to measure convolutions of GPDs, lattice QCD allows us to extract hadron 
matrix elements for the twist-2 operators, which can be expressed in terms of generalized form factors. 

In order to compare hadron matrix elements of these local operators to experiment one needs to renormalize them. 
The aim of this paper is to calculate non-perturbatively the renormalization factors of the above twist-two fermion 
operators within the twisted mass formulation. We show that, although the lattice spacings considered in this work are 
smaller than 1 fm, O(a^) terms are non-negligible and significantly larger than statistical errors. We therefore compute 
the 0(a^)-terms perturbatively and subtract them from the non-perturbative results. This subtraction suppresses 
lattice artifacts considerably depending on the operator under study and leads to a more accurate determination of 
the renormalization constants. Preliminary results of this work have been published in Refs. [11, 12]. 

The paper is organized as follows: in Section II we give the expressions for the fermion and gluon actions we 
employed, and define the twist-two operators. Section III concentrates on the pcrturbative procedure, and the 0(0^)- 
corrected expressions for the renormalization constants Zq and Zo- Section IV focuses on the non-perturbative 
computation, where we explain the different steps of the calculation. Moreover, we provide the renormalization 
prescription of the RI'-MOM scheme, and we discuss alternative ways for its application. The main results of this 
work are presented in Section V: the reader can find numerical values for the Z-factors of the twist-2 operators, which 
are computed non-perturbatively and corrected using the perturbative 0(0^) terms presented in Section III. Since in 
general Z-factors depend on the renormalization scale, we also provide results in the RI'-MOM scheme at a reference 
scale, fi ^ 1/a. For comparison with phenomenological and experimental results, we convert the Z-factors to the MS 
scheme at 2 GeV. In Section VI we give our conclusions. 

A forthcoming paper [13] will focus on the perturbative procedure and will present results for local fermion operators 
(scalar, pseudoscalar, vector, axial, tensor). 
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II. FORMULATION 



Lattice action 



For the gauge fields we use the tree- level Symanzik improved gauge action [14], which includes besides the plaquette 
term J7^_^^y also rectangular (1 x 2) Wilson loops U]^^^^ 



^P = fE(^o E ^ {l-ReTr(C/ifJ} 



(3) 



l<p<i/ 



with /3 = 2 Nc/g^, bi = —1/12 and the (proper) normalization condition bg — 1 — 8bi. Note that at 6i = this action 
becomes the usual Wilson plaquette gauge action. 

The fermionic action for two degenerate flavors of quarks in twisted mass QCD is given by 



Sf = a'^'^x{x){Dw[U] + mo + iHol5T^)x{x) 



(4) 



with the Pauli matrix acting in the isospin space, /xq the bare twisted mass and Dw the massless Wilson-Dirac 
operator defined as 



where 



^Mx) = - 



U^{x)il){x + aft) - il){x) 



and '^*tp{x) 



UUx — afj,)tp{x — ajl) — tp{x) 



For completeness we also provide the definition of the backward derivatives 



il){x)^ n = - il){x + afi)U^{x) — ■tp{x) and tl){x)^*^ = — - 



tjj{x — afb)Ufj,{x — ajl) — tlj{x] 



(5) 



(6) 



(7) 



Maximally twisted Wilson quarks are obtained by setting the untwisted bare quark mass toq to its critical value mcr, 
while the twisted quark mass parameter /io is kept non-vanishing in order to give the light qiiarks their mass. In 
Eq. (4) the quark fields x are in the so-called "twisted basis" . The "physical basis" is obtained for maximal twist by 
the simple transformation 



tp{x) = exp ( j x{x), tp{x) = x{x) exp ( ^Isr^ 



(8) 



In terms of the physical fields the action is given by 

St = E^(^) + - n^r' (-y + mc.) + Mo) ^(x) . 

In this work we consider twist-two operators with one derivative, which are given in the twisted basis as follows 



(9) 



{nv} _ 



DV 



= X7{/i D ^}T«x 



= < 



(10) 



o 



DA 



.V'757{m-D,.}t3V' 



a = 1 
a = 2 
a = 3 



(11) 



O 



tJ.{v p} 
DT 



Xl5(^t,{-^ D pyT°-X 



^75CT^{i. ^}t"V a = 1,2 
-iijja^{^ D pylip a = 3 



(12) 
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with the covariant derivative defined as 

s=i[(%ly-%lj)]. (13) 

The above operators are symmetrized over two Lorentz indices and are made traceless 

oi'^^y = -(0"^ + 0'"'^ --S'^'^^o^^ . (14) 

A 

This definition avoids mixing with lower dimension operators. We denote the corresponding Z-factors by Z^y, Z^^^, 
Z^rj, . In a massless ronormalization scheme the rcnormaUzation constants arc defined in the chiral fimit, where isospin 
symmetry is exact. Hence, the same value for Z is obtained independently of the value of the isospin index a and 
therefore we drop the a index on the Z-factors from here on. However, one must note that, for instance, the physical 

V'7{ju^i/}T^V is renormalized with Z-qa, while '^7{^^^}r^V requires the .^dVj which differ from each other even in 
the chiral limit. 

The one derivative operators fall into different irreducible representations of the hypercubic group, depending on 
the choice of indices. Hence, we distinguish between 

(15) 
(16) 

(17) 
(18) 
(19) 
(20) 

Thus, Zuvi will be different from Zuv2, but renormalized matrix elements of the two corresponding operators will be 
components of the same tensor in the continuum limit. Odti is sufficient to extract all generalized form factors of the 
one derivative tensor operators. Although, in this work, we will calculate non-perturbatively only the renormalization 
constants for the vector and axial- vector operator, we will provide the perturbative O(a^) terms also for the tensor 
operators. 





= Odv 


with 


/« = 


v 






= Odv 


with 




v 




Odai 


= 0DA 


with 


M = 


V 




CdA2 


= OuA 


with 




V 




Cdti 


= Out 


with 




V 


= p 


ObT2 


= Out 


with 




V 





III. PERTURBATIVE PROCEDURE 



Here we present our results for the renormalization factor of the fermion field, Zq, that enters the evaluation of the 
twist-2 operators, and we also calculate Zu\r, Zua, Zut- Our calculation is performed in 1-loop perturbation theory 
to 0(0,^). Extending the calculation up to ©(a^) brings in new difficulties, compared to lowers order in a; for instance, 
there appear new types of singularities. The procedure to address this issue is extensively described in Ref. [15]. Many 
IR singularities encountered at C(a^) would persist even up to 6 dimensions, making their extraction more delicate. 
In addition to that, there appear Lorentz non-invariant contributions in 0{a^) terms, such as J2fj.Pti/p'^ ^ where p is 
the external momentum; as a consequence, the Z-factors also depend on such terms. The knowledge of the order 
terms is a big advantage for non-perturbative estimates, since they can eliminate possible large lattice artifacts, once 
the O(a^) perturbative terms are subtracted. 

For all our perturbative results we employ a general fermion action, which includes the clover parameter, csw, and 
non-zero Lagrangian mass, m. For the fermion field renormalization we also have a finite twisted mass parameter, /xQ) 
so we can explore its /xo dependence. Only the renormalization factors of the twist-2 operators are obtained at /xq = 0, 
but we still consider m =^ 0. For gluons we use Symanzik improved actions (Plaquette, Tree-level Symanzik, Iwasaki, 
TILW, DBW2) [15]. The purpose of using such general fermion and gluon actions is to make our results applicable to 
a variety of actions used nowadays in simulations. The expressions for the matrix elements and the Z-factors are given 
in a general covariant gauge, and their dependence on the coupling constant, the external momentum, the masses and 
the clover parameter is shown explicitly. 
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A. Fermion field renormalization 



The fermion action used for the inverse fermion propagator improvement is a combination of Wilson/clover/twisted 
mass fermions, given by 



(21) 



All quantities appearing in Eq. (21) are defined in Ref. [15]. There are two 1-loop diagrams involved in this particular 
computation, which are illustrated in Fig. 1. 




1 2 

FIG. 1: One-loop diagrams contributing to the fermion propagator. Wavy lines represent gluons and solid lines fermions. 
We compute Zq in the RI'-MOM renormalization scheme, defined as 



1. 



-Tr 



-Tr 



Pp=fj-P 



|Ep7pSin(aPp) 



(22) 



Pp=fj-p 



For comparison reasons, an additional renormalization prescription was also applied (see Eqs. (56) - (57)). Since we 
want to take into account all 0{a^) terms, we perform a Taylor expansion leading to 



Ep7p(Pp 



6 i^pJ 



T.pPl 



3 EpPl 



' '5'l-loop(p) 



+ 0(a'g^g') 



Pp=Pp 



I 

-Tr 

4 



^2 ' '^l-loop(P) Q o ^2 



1 / 



2p2 (^2)2 



'5'l-loop(p) 



(23) 



Pp=Pp 



The trace is taken only over spin indices and '5']^_3jqqp is the inverse fermion propagator that we computed up to 1-loop 
and up to 0{a^). We make the following definitions for convenience: = J2pP% p4 = J^pPp' ^ — J2plpPp 

= ^plpP^- A very important observation is that the 0{a^) terms depend not only on the magnitude, X^^Ppi 
but also on the direction of the external momentum, pp, as manifested by the presence of the terms ^pPj,- As 
a consequence, alternative renormalization prescriptions, involving different directions of the renormalization scale 
Mp — Pp, treat lattice artifacts diversely. 

For the special choices: csw = 0, r = 1 (Wilson parameter), A = (Landau gauge), mo — 0, fiQ — 0, and for 
tree-level Symanzik gluons, Zq can be read from the following trace 



1 

4^^ 



StrLiP) ■S'lAoopb) 



Pp=p.p 



l+g^|-13. 02327272(7) 



73 



m2(i.14716212(5)-— In(aV')) 



4 

pf^P , 



1 ^7 n ~i 

2.1064995(2) - — \n{a' M^))] } + 0{a^ g\g') (24) 



where g^ = g^ C p / {IQt^^) , Cp ^ [N^ - l)/(2iVc) and /^^ ^ Y^p^-l- In Ref. [13] we provide Zq for csw = 0, A = 
and the dependence on mg, /^o shown explicitly. Its most general expression (including csw and A, as well as a wider 
choice of values for the parameters entering the Symanzik action) is far too lengthy to be included in paper form (Zg, 
Zdvj ^da, ^dt around: 250, 800, 800, 950 terms respectively); it is provided in electronic form along with Ref. [13]. 
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B. Renormalization of twist-2 operators 



Here we present the computation of the amputated Green's functions for the following three twist-2 operators 



O 



DV 



o 



DA 



O- 



DT 



T 



(25) 
(26) 
(27) 



which, being symmetrized and traceless, have no mixing with lower dimension operators. 

The Feynman diagrams that enter our calculation are given below, where the insertion of the twist-2 operator is 
represented by a cross. We have computed, to O(a^), the forward matrix elements of these operators for general 




FIG. 2: One-loop diagrams contributing to the computation of the twist-2 operators. A wavy (solid) line represents gluons 
(fermions). A cross denotes an insertion of the operator under study. 



external indices /i, v (and p for the tensor operator), external momentum p, toq, g, N^, a, csw ^^'^ ^- Setting A = 1 
corresponds to the Feynman gauge, whereas A = corresponds to the Landau gauge. Our final results were obtained 
for the 10 sets of Symanzik coefficients we have used in the calculation of Zq [15]. 

In order to define Zq, we have used a renormalization prescription which is most amenable to non-perturbative 
treatment: 



Zr^ZoTr 



Tr 



(28) 



where L'-' denotes the amputated 2-point Green's function of the operators in Eqs. (25) - (27), up to 1-loop and up 
to 0(0^). The tree-level expressions of the operators including the 0{a^) terms are 



7m - « 



.2 Pi- 



2 



+ 0(a^) 



.2 P. 



-DA2 



(p)=75ij^rb) 



(29) 

(30) 
(31) 



rDT2^ \ 
-t^tree VP) 



2 75 



<y,u (pp - f ) + ^PP (p- - f ) ) + ^i^') 



(32) 



We perform a Taylor expansion up to 0{a?) in the right hand side of the renormalization condition and it leads to 
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the following 



Tr L^^IHp) ■ KL\P) 



Tr 



Tr 



12 



12 



(33) 
(34) 

(35) 



Tr 



-DT2/ 



(p) • Lll^ip) =pI+pI-'1^ (pI+p^^) + 0{a') 



(36) 



For the special choices: mo = 0, csw = 0, r = 1, A = (Landau gauge), and for tree-level Symanzik gluons, we 
obtain for the left hand side of Eq. (28) 



Tr[L°-^(p).LSr(p)] 



+ f{\^ +p' (3.610062(3) - |ln(a^p^)) + pi (27.54716(3) - ^ln(a^p2)) 



[(p^)^ (0.11838(2) + ^ln(a>^)) +/p^ (-0.6573(1) - ^ln(a^p^)) 
+ p4 (-1.71886(3) + |Iln(.^pV^§) 



^ -1/ i«in.o^t;N^94 2 2 29 p4 ^ 169 p' ] T 
+ p, (-16.1049(5) + -ln(a p ) + - ^ + __)J | 



90 (p2 



45 p2 - 



(37) 



Tr[L°^^(p).LSr(p) 



+ 



2 2 ^ / 4 4\ 

-Pm ^P'^ + "3"(Pm 

^1^4^ + (p' +P') (15.04575(1) - |ln(a^/)) 
+ a' [(p^+pf.) (-7.1429(1) + |iln(a^p^)) 

+ (pI+pD (/ (-0.13212(3) - 1^ ln(.r/)) + g^^) 



_L 2 2. ,nnaRn^_Ll013 2 2.^29 p4 , 169 (p^ + p^) ^] ] 
+ p,p. ( - 4.0096(1) + — Ha p) + ^J-^ + ^^^)J } 



90 (p2 



90 



+ 0(a',/) 



Tr[L°^^(p).L°i^(p)] = 2p^ + i/ + a^(-lp4-|pt) 



+ 3'{-^4 (-4.127332(3) + |ln(a-p^)) +p^ (-31.68532(3) + ^ln(a^p^)) 



i p- 
+ a 



3 

[{pY (0.17035(2) + ^ln(a>^)) +/p^ (0.3982(1) - ^ln(a^p^)) 



288 
397 
720 
2 

5 



+ p4(l.69230(3)-||iln(a^/) + ^|) 

■ 4/1c/,«1Q/c^ , 2 2^ 29 p4 169pL]\ 

+ p, (18.4613(5) + -ln(a p ) - - ^ - __)J| 



+ 0(a',/) 



(38) 



(39) 
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It \ T DA2 / nI 2,2 4 

(P) ■ -t'troc {P)\ = Pn+Pv- y (Pm 



+ (^'^ ( - 16.10196(1) + |ln(a^p^)) 

491 



+ a' [(Pm + Pt) (7.2286(1) - ln(a^ p^)) 

+ (P^ +P?.) (/ (0.75869(3) - i|ln(a^/)) - |^^) 

_L 2 2 ^ 187, , 2 2. 29 p4 169 (p^fp^.]-, 

+ p,p. (4.8509(1) + — ln(a p ) - - ^ - ^^^)J } 

+ 0{a\g*) (40) 



rp frDTl/ X rDTl/ n1 P P/j , n 2 , 2 /P/J P^ 2p* 
Tr \L (p) ■ it.ee (P)J = ^ + 2p^ + (I (-^^-j^ —) 



+ f [ipl-p') (4.226559(3) - \n{a p^)) +pl{- 29.11666(2) + 5 In(a^p^)) 

43 



+ a^[p^(- 0.14754(2) -^ln(a^p^)) 



+ p4 (1.93789(3) - — ln(a p )-__ + _ _) 



+ p' (l.7215(l)p?. + ^ liV p^)p^ +0.37022(2) p?, - ^ ln(a^ p2)p^) 

+ p^ (14.9155(4)- 11 ln(a^/)-^|) 

+ p?p?. (2.4896(1) -|iln(.r/)- 

+ pt(- 2.24911(4) + Ziln(a^/))-^0]} 
+ 0(a',/) (41) 



4 4 

rrn f r DT2 / N r DT2 / nI 2 , 2 , 2 / Pp Pi/ X 

Trp (p) • itrco (P)j = Pp+P,. + a (~3 



1 -2 

+ g 



{{pI+pI) ( - 15.84740(1) + 3 ln(a2p")) 

+ a' [{pi +pI) (0.22134(3) p' + ^ ln(a^p^)p^ - ^ ^ 



+ O.736O4(2)p?.-|iln(.r/)p?.-^0^ 



67 PpP^P^ , / 4 , 4\ on/in/i N 1051 , . 2 2n\ 

+ (P.+Pp) (7.3949(1) - — — ln(a p )) 



15 p2 ■ ■ ^ 720 

1609 2 __2^^ 67 PpP^_+PppJl 

p2 

+ 0{a\g^) (42) 



, 2 2/ono.t-n/o^ 1609 1 t 2 2xx 67 PpP^+PpP^ -\\ 
+ PpP, (2.98450(8) - — ln(a P )) " ^ -5 J | 



In our forthcoming publication [13], wc will include an ASCII file with all our numerical results for general value of 
^7 cswj ™o a-nd the 10 sets of Symanzik gluon actions; the file is best perused as Mathematica input. In addition, 
there appears an Appendix providing the exact 0{a^) terms that need to be subtracted from Zq (a fraction of the 
two traces and Zq of Eq. (28)). 

The 0(0^) terms shown in Eqs. (37) - (40) arc used to correct our non-pcrturbative results for Zdvi, ^dv2! Zj^ai, 
Zda2, in order to better control artifacts. For the subtraction procedure we use the boosted coupling [16] instead 
of the bare one 



^boosted ~ /„, \ ' (4^^) 
y^plaq ) 



where (upiaq) is the plaqucttc mean value. 
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IV. NON-PERTURBATIVE CALCULATION 

A. Evaluation of correlators 

In the literature there are two main approaches that have been employed for the non-perturbative evaluation of the 
renormalization constants. They both start by considering that the operators can all be written in the form 

0{z)=Y,u{z)J{z,z')d{z'), (44) 

z' 

where u and d denote quark fields in the physical basis and J denotes the operator we are interested in, e.g. J{z^ z') = 
^z,z'liJi would correspond to the local vector current. For each operator we define a bare vertex function given by 

„12 _ 

^(^') = ir E e-'''^''-''\u{x)u{z)J{z,z')d{z')d{y)), (45) 

x,y,z,z' 

where p is a momentum allowed by the boundary conditions, V is the lattice volume, and the gauge average is 
performed over gauge-fixed configurations. We have suppressed the Dirac and color indices of G{p). The first 
approach relies on translation invariance to shift the coordinates of the correlators in Eq. (45) to position z = Q 
[17, 18]. Having shifted to z = allows one to calculate the amputated vertex function for a given operator J for any 
momentum with one inversion per quark flavor. 

In this work we explore the second approach, introduced in Ref. [19], which uses directly Eq. (45) without employing 
translation invariance. One must now use a source that is momentum dependent but can couple to any operator. 
For twisted mass fermions, we use the symmetry S'^{x,y) = ^58'^^ {y, x)^5 between the u— and d— quark propagators. 
Therefore with a single inversion one can extract the vertex function for a single momentum. The advantage of this 
approach is a high statistical accuracy and the evaluation of the vertex for any operator including extended operators 
at no significant additional computational cost. Since we are interested in a number of operators with their associated 
renormalization constants we use the second approach. We fix to Landau gauge using a stochastic over-relaxation 
algorithm [20] , converging to a gauge transformation which minimizes the functional 

F = 5^ Re tr [U^{x) + U^x - A)] • (46) 

Questions related to the Gribov ambiguity will not be addressed in this work. The propagator in momentum space, 
in the physical basis, is defined by 

= ^ E H^My)) , s'^ip) = ^Y1 ^"^"^'"'^ {d{x)d{y)) ■ (47) 

x,y x,y 

An amputated vertex function is given by 

Tip) = iS-{p))-'G{p){S\p)r'. (48) 
and the corresponding renormalized quantities are 

Snip) = ZgS{p) , Tnip) = Z-^ZoT{p) , (49) 

In the twisted basis at maximal twist, Eq. (45) takes the form 

^(P) = 1^7 E e-'f (^-«) ((1 + i^5)u{x)u{z){t + 115) J{z, z'){t - ij5)d{z') d{y){l - ns)) . (50) 

x,y,z,z' 

After integration over the fermion fields, and using S^{x,z) = 75>S''\^;, a;)75 this becomes 

„12 / ^ 4. \G 

Gip) = ^J^(iT^-^l5)S^ iz,p){t-ij,)Jiz,z'){l-ij,)S^iz,p){t-tj,)j , (51) 

z ^ ' 

where (...)'^ is the integration over gluon fields, and S{z,p) = ^ye.'^^^S{z,y) is the Fourier transformed propagator 
on one of its argument on a particular gauge background. It can be obtained by inversion using the Fourier source 

bl{x) = e'^-5^05ab, (52) 
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for all Dirac a and color a indices. The propagators in the physical basis given in Eq. (47) can be obtained from 

Z 

^"(P) = -jEe+^''((l-*^5)5'^\^,p)(l-*75))^, (53) 

z 

which evidently only need 12 inversions despite the occurrence of both u and d quarks in the original expression. 

We evaluate Eq. (50) and Eq. (53) for each momentum separately employing Fourier sources over a range of a^p^ 
for which perturbative results can be trusted and finite a corrections are reasonably small. 

B. Renormalization Condition 

The renormahzation constants are computed both perturbatively and non-perturbatively in the Rl'-MOM scheme at 
different renormalization scales. We translate them to the MS-scheme at (2 GeV)^ using a conversion factor computed 
in perturbation theory to 0{g^) as described in Section V. The Z-factors are determined by imposing the following 
conditions: 



(5^(p))-iS(°)(p)l ^ ^ (54) 



r^.(p)r(S-^(p)l = 1, (55) 



where is the renormalization scale, while Sl and Fj, correspond to the perturbative or non-perturbative results. 

The trace is now taken over spin and color indices. These conditions are imposed in the massless theory, i.e. at 
critical mass and vanishing twisted mass. At finite lattice spacing there are two choices for 5*^°-' and F'''' entering 
Eq. (55). One can take either the tree level or the continuum results for 5^°^ and F(°\ which differ by 0(a^)-terms. 
The continuum free propagator in terms of continuum momentum is 

5(0)(p) = (56) 

fW(p) = -iOi^p^y. (57) 

We refer to this choice as method 1. A different choice is to define the free propagator using the lattice momentum [18, 
19] : 

5<o.,rt = (58) 
EpSm(pp)2 

fW(p) = -id{^ sin(p,}), (59) 

which we will refer to as method 2 used in Ref. [21]. 

The choices for S^'^^ and F^*]] given in Eqs. (58) - (59) arc preferable, compared to Eqs. (56) - (57), since only for 
method 2 we obtain Zq = l, Zq = 1 when the gauge field is set to unity. Similarly, in the perturbative computation 
only method 2 gives Zq = \, and Zq = 1 at tree-level. On the contrary, the Z-factors obtained from method 1 have 
lattice artifacts even at tree- level. In this sense, method 2 is an improvement of method 1. and can be thus considered 
superior. Obviously, the renormalization constants using the two methods differ only in their lattice artifacts, as 
can be seen by Eq. (23). We find that, for the cases considered here, non-perturbative results using method 2 lead 
to Z-factors with smaller lattice effects. We demonstrate this by examining the following case: Let us consider the 
momenta as given in Table 11, which fall into two sets, those with spatial components 2 tt (2, 2, 2)/i and those with 
2 TT (3, 3, 3)/L; there is only one non democratic momentum, with spatial components 2 tt (3, 3, 2) / L, but this behaves 
similarly to the second set mentioned above. The two sets of momenta do not fall on the same curve, a behavior that 
is due to cut-off effects. This is clearly seen in Fig. 3 where we present Zj^p^ at /3 = 3.9 and /txo = 0.0085 using the 
two methods (upper plot) (similarly for Z^y in Fig. 4). The statistical errors in Fig. 3 and the rest of the graphs are 
smaller than the size of the symbols. As can be seen, the two sets of momenta differ in particular when using method 
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FIG. 3: Z^''^^ for 13 = 3.9 (a"i=2.217 GeV) and = 0.430 GeV for method 1 (open symbols) and method 2 (filled symbols). 
The upper plot corresponds to non-perturbative results, where the index A, B represents the set of momenta with spatial 
components 2 tt/L (3, 3, 3) and 2 tt/L (2, 2, 2), respectively. The lower plot shows the non-perturbative results after subtracting 
the perturbative 0{g^ a^)-terms, where the two methods give almost identical results. Moreover, in method 1, the jump between 
the two sets of momenta disappears. 





FIG. 4: In the left panel we show Zf^ using the same notation as in Fig. 3. In the right panel, upper graph we show Z^^ 
again using the notation of Fig. 3 whereas in the lower graph we show a comparison between method 1 after subtracting the 
perturbative C(a^)-terms (diamonds) and method 2 without any subtractions (filled squares). 



1. However, it is important to note that, after subtracting the 0{a^) perturbative contributions one obtains values 
that are consistent between the two methods (see lower plot). Moreover, in method 1 the jump observed between 
the two sets of momenta disappears. For method 1 the subtraction of the perturbative O(a^) terms refers to both 
contributions of 0{g^ a^) (tree-level) and 0{g^ a^) (1-loop). Thus, upon subtraction method 1 can be compared to 
method 2 because we remove all C(a^) terms up to 1-loop, leading to almost equivalent results as can be seen in the 
lower plot in Fig. 3. We would like to stress that the substraction is necessary and yields results superior to using 
the unsubtracted method 2. This is demonstrated in Fig. 4 in the right hand size plots where in the upper plot we 
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show the unsubtracted methods 1 and 2, while in the lower plot we show a comparison method 1 after subtraction of 
perturbative 0(a^)-terms with the unsubtracted method 2. One can observe that the jump between different set of 
momenta, appearing in the unsubtracted mctliod 2, almost disappears in the subtracted method 1. The same pattern 
appears in all Z-factors of the one-derivative operators, as well as for Zq. The latter has an impact on all other 
renormalization constants discussed here. This effect, as expected, becomes less pronounced at /3 = 4.05 and 4.20, 
and disappears for small a^p^ as demonstrated in the next section. The results presented in all Tables correspond to 
the Z-factors obtained using method 2. 

V. RESULTS 
A. RI'-MOM condition 

We perform the calculation of renormalization constants for three values of the lattice spacing corresponding to 
/3 = 3.9, 4.05 and 4.20. The lattice spacing as determined from the nucleon mass is 0.089 fm, 0.070 fm and 0.056 fm 
respectively. For ,5 = 3.9 we consider three different quark masses, corresponding to = 0.302 GeV (a/io = 0.004), 
= 0.376 GeV (a/xo = 0.0064) and m.^ = 0.430 GeV (a/Uo = 0.0085), in order to explore the dependence of the 
Z-factors on the pion mass. At ^ = 4.05 we consider two volumes, 24^ x 48 and 32^ x 64 in order to check for 
finite volume effects. To extract the renormalization constants reliably one needs to consider momenta in the range 
^QCD < p < 1/a. We relax the upper bound to be ~ 2/ a to 3/a, which is justified by the linear dependence of 
our results on a^. Therefore, for each value of /3 we consider momenta spanning the range 1 < a^p^ < 2.7 for which 
perturbation theory is trustworthy and lattice artifacts are still small enough. In Table I we summarize the various 
parameters of the action, that we used in our simulations. 





a (fm) 




(GoV) 




X T 


3.9 


0.089 


0.0040 


0.3021(14) 


24^ 


X 48 


3.9 


0.089 


0.0064 0.37553(80) 24^ 


X 48 


3.9 


0.089 


0.0085 


0.4302(11) 


24'^ 


X 48 


4.05 


0.070 


0.006 


0.4082(31) 


24'^ 


X 48 


4.05 


0.070 


0.006 


0.404(2) 


32'^ 


X 64 


4.05 


0.070 


0.008 


0.465(1) 


32^ 


X 64 


4.20 


0.055 


0.0065 


0.476(2) 


32^ 


X 64 



TABLE L Action parameters used in the simulations. 

In Table II we present the statistical sample for the parameters and momenta we used in the simulations. Using the 
number of configurations shown in Table II leads to results with very high statistical accuracy, easily below 0.5%. 
The results for the subtracted Z-factors (method 2) at /3 = 3.9 are tabulated in Table III for the highest and lowest 
twisted mass parameter used (for the lowest mass we have obtained the Z-factors only for 6 momenta). Comparison 
between the Z— factors for two different masses shows that any dependence on the pion mass is within the small 
statistical errors. This negligible dependence is not a result of the O(a^) subtraction, as demonstrated in Fig. 5. 
The left plot illustrates the pion mass dependence of the unsubtracted .^dvi for three renormalization scales ranging 
from 5.75 GeV^ to 11.75 GeV^, while the subtracted ^dvi is shown in the right plot. The same behavior is observed 
for all renormalization constants considered here. The subtracted Z-factors (method 2) for /3 = 4.05 and /3 = 4.20 
are presented in Tables IV-V, respectively. In order to see possible volume effects we compute the renormalization 
constants at /3 = 4.05, fio = 0.006, for two lattices with different size, namely for 24"^ x 48 and for 32^ x 64. For 
this comparison we used momenta that correspond to the same renormalization scale: For the small lattice we use 
27r(3/48, 3/24, 3/24, 3/24), in lattice units, whereas for the larger one we employ 27r(4/64, 4/32, 4/32, 4/32). The 
volume effects appear to be ~ 0.1%, as can be seen from Table VI. 



Given the small statistical errors one may carefully examine the systematic errors. As already noted a systematic 
effect comes from the choice of S^°^ and r^°\ To give an example, at /3 = 3.9, /Uq = 0.004, /j,'^ w 5.75 GeV^ method 1 
leads to Zq = 0.76606(7) while method 2 gives Zq = 0.80514(7), before any subtraction of 0(0^) is carried out. This 
systematic effect is removed after perturbative subtraction is applied. Another, much smaller, systematic effect comes 
from the asymmetry of our lattices both because they are larger in their time extent and because of the antiperiodic 
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TABLE II: Statistical sample at /? = 3.9, 4.05, 4.20 for various momenta. 
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FIG. 5: Zovi at /3 = 3.9, as a function of the pion mass: = 0.302 GeV (a/to = 0.004), = 0.375 GeV {a/io = 0.0064) and 
m,r = 0.429 GeV (a/io ~ 0.0085). The left plot regards the unsubtracted non-perturbative results and the right one corresponds 
to the subtracted data. 



boundary conditions in the time direction. For instance, using the same /3 and /io as in the previous example, method 
1 at W7.6 GeV'^, in the temporal direction of the current gives Zdvi — 1.1387(2) while the average from the 
three spatial directions leads to .^dvi = 1-1006(2). This effect can be seen in Fig. 6 where we plot separately the 
renormalization constant Zdvi determined from the temporal indices, the spatial indices and the average of those 
two. In the same figure we also show that upon subtraction this systematic effect disappears (lower plot). For 
Tables III - V we use for Zuvi the average of Z^yy, Z^y with j/ = 1,2,3, while for Zdv2 the average of Z^^y, Z^y 
with 7^ p = 1, 2, 3. We apply the same procedure for the twist-2 axial operator. 

Chiral extrapolations are necessary to obtain the renormalization factors in the chiral limit. As already pointed out 
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TABLE III: The renormalization constants at ^ = 3.9 with /Uo = 0.004, 0.0085 for lattice size: 24^ x 48. 
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TABLE IV: Renormahzation constants at /3 = 4.05, afxo = 0.008 for lattice size 32^ x 64. 

the dependence on the pion mass is insignificant. Allowing a slope and performing a linear extrapolation to the data 
shown in Fig. 5 yields a slope consistent with zero. This behavior is also observed at the other /3-valucs and therefore 
the renormalization constants are computed at one quark mass, given in the Tables III-V. Figures 7, 8, 9 demonstrate 
the eflFect of subtraction, for all three ^ values, as a function of the renormalization scale (in lattice units). For all 
cases we observe a significant correction upon subtraction; the lattice artifacts for ^da2 turn out to be very small 
for most values of the momentum. In addition, the lattice artifacts decrease by employing higher values for /3 (finer 
lattice), as expected. 

B. RI'-MOM at a reference scale 

All our Z-factors have been evaluated for a range of renormalization scales. In this subsection we use 2-loop 
perturbative expressions to extrapolate to a scale /x = I/o (the values for a are taken from Table I). Thus, each result 
is extrapolated to 1/a, maintaining the information of the initial renormalization scale at which it was compTitcd. 
Although the 3-loop formula is available for the following expressions, the 0{g^) corrections are insignificant compared 
to the lower order results. 

The scale dependence is predicted by the renormalization group (at fixed bare parameters), that is 

Zgi'(M) = i?o(M,Mo) Zl^\iia) (60) 
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FIG. 6: (squares), Z^-^ (circles), Z'g^^ (crosses), for ^ = 3.9 (a-^=2.217 GeV), 



0.430 GeV using method 1. The 



upper plot corresponds to the purely non-perturbative results, while the lower plot shows the non-perturbative results after 
subtracting the perturbative terms of C(a^). 




FIG. 7: Renormalization scale dependence for the Z-factors at P — 3.9 and m-^ — 0.430 GeV 
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FIG. 8: Renormalization scale dependence for the Z-factors at /3 = 4.05 and — 0.465 GeV 
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FIG. 9: Renormalization scale dependence for the Z-factors at /3 = 4.20 and — 0.476 GeV 



17 



[nt ,nx,ny,nz) 


^DVl 




.^DAl 




(4,2,2,2) 


1.1585(4) 


1.215(1) 


1.2266(5) 


1.257(1) 


(5,2,2,2) 


1.1387(4) 


1.189(1) 


1.2052(5) 


1.230(1) 


(6,2,2,2) 


1.1203(3) 


1.1642(9) 


1.1853(4) 


1.2040(9) 


(3,3,3,2) 


1.1069(2) 


1.1459(9) 


1.1702(3) 


1.1855(9) 


(7,2,2,2) 


1.1028(2) 


1.1413(8) 


1.1668(3) 


1.1804(8) 


(2,3,3,3) 


1.0943(2) 


1.1257(8) 


1.1599(3) 


1.1643(8) 


(8,2,2,2) 


1.0853(1) 


1.1197(4) 


1.1488(2) 


1.1579(4) 


(3,3,3,3) 


1.0841(1) 


1.1177(8) 


1.1438(3) 


1.1577(8) 


(4,3,3,3) 


1.0743(2) 


1.1079(7) 


1.1293(2) 


1.1491(7) 


(5,3,3,3) 


1.0651(2) 


1.0968(6) 


1.1163(2) 


1.1390(6) 


(6,3,3,3) 


1.0511(2) 


1.0812(5) 


1.0983(3) 


1.1274(5) 


(10,2,2,2) 


1.0528(1) 


1.0856(3) 


1.1189(1) 


1.1223(3) 



TABLE V: Renormalization constants at /3 = 4.20, /io = 0.0065 for lattice size: 32^ X 64. 



lattice Zuvi Zr>v2 Zuai Z0A2 

24^x48 1.0700(2) 1.0923(2) 1.1190(2) 1.1117(2) 
32^x64 1.07123(6) 1.0928(2) 1.12037(7) 1.1122(2) 



TABLE VI: Renormalization constants at /3 = 4.05, /Jo = 0.008 using method 2 and two lattice sizes: 32^ x 64 for (4,4,4,4) and 
24^x48 for the rest of the momenta. 
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To 2 loops, the running coupling, /3-function and anomalous dimension 7 are as follows: 
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The expressions for the anomalous dimension of the fermion field and the twist-2 vector/axial operators are given 
in Ref. [22], 
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FIG. 10: Renormalization factors in the RI'-MOM scheme at renormalization scale 1/a, for /3 = 3.9, /io — 0.0085. The black 
circles correspond to the unsubtracted results, while the magenta diamonds to the results with perturbatively subtracted one 
loop O(a^) artifacts. 



where Tp = 1/2, Ca = N^. Using Eqs. (60) - (66) we obtain the Z-factors at ^ = 1/a for /3 = 3.9, 4.05, and 4.20, 
which are plotted in Figs. 10 - 12. 



C. Conversion to MS 



The passage to the continuum MS-scheme is accomphshed through use of a conversion factor which is computed 
up to 3 loops in perturbation theory. By definition, this conversion factor is the same for the vector and axial twist-2 
renormalization constant, but will differ for the cases Zdvi (^dai) and Zdv2 (2^dA2)i that is 



Cdvi = Cdai = pjY (67) 

■^DVl 

Cdv2 = CdA2 = j^p ■ (68) 

'^DV2 

This requirement for different conversion factors results from the fact that the Z-factors in the continuum MS-scheme 
do not depend on the external indices, /x, (see Eq. (2.5) of Ref. [22]), while the results in the RI'-MOM scheme do 
depend on fi and h'. Of course the conversion factors take a different value for each renormalization scale; actually, 
the direction of the momentum is required to be known (Eqs. (85) - (86)). 

The 3-loop expressions for the conversion factors from our RI'-MOM scheme (Eq. (55)) to the MS do not appear 
directly in the literature, but can be extracted using results from Ref. [22]. In the latter publication the reader can find 
the conversion factor from an alternative definition of RI'-MOM (which we denote by RI) to the usual MS, G^^^j~iu^- 
This alternative definition reads 
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FIG. 11: Same as Fig 10, but for /3 = 4.05 and fio = 0.008. 



where T,'^^'> can be extracted from the bare amputated Green's function as follows 



(70) 



The author of Ref. [22] provides the 3-loop expression for the renormalized S*^^) in the scheme of Eq. (69) (note that 
by definition the renormalized E'"'^' equals 1 at = fi^). These elements can be used to reconstruct the renormalized 
Green's function 



(2) RI' finite, 



(p)^ 



(71) 



in which we apply our RI'-MOM condition in order to obtain 



Z— ^ ^Ri 



yRl' 



(72) 
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FIG. 12: Same as Fig 10, but for P = 4.20 and = 0.0065. 



Once we have these two elements we extract the conversion factor of Eqs. (67) - (68) up to 3 loops, 



CDVi(Ai) = 

CdV2(m) = 



^RI 

'^DVl 
^RI 

7RI' 



^RI 

^rF 

^DVl 

^ri 
7RI' 

'^DV2 



^DV 



^RI 



^RI 



(73) 
(74) 



The conversion to the MS is then given by 



and 



7MS 

^DV2 

MS 
^DA2 



in) = Cdv2(m) ■ 2'5v2(Ai) 
Z^^2{^^) = C'dV2(m) ■ -^1)12 (m) : 



(75) 
(76) 

(77) 
(78) 



which correspond to the Z-factors at the same renormalization scale in the RI'. One wants to obtain the renormal- 
ization constants at the scale of 2 GeV, and to do this we use the 2-loop formula in Eq. (61)-(62) to evolve the scale 
from /i to 2 GeV. In these formulas we need to insert the anomalous dimension in the MS-scheme which read [22] 



^Ja) = 2 ^Cpa + 2 ^ [47C7^ - UC'f - mTpMp] a' 



(79) 



where a — j (167r^). The additional factor of 2 that we included, comes from the different definition of the anomalous 
dimension that leads to Ref. [23]. To summarize, the Z-factors in the continuum MS -scheme at = 2 GeV are given 
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For the SU{Nc — 3) colour group {Ca = 3, Cp — 4/3, Tp = 1/2), Landau gauge (A = 0), and general quark 
flavours, we have the following conversion factors 
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where a = f;^/(167r^) and C(") is the Riemann Zeta function. 

A "renormalization window" should exist for Ag,^^ << /i^ << 1/a^ where perturbation theory holds and finite-a 
artifacts are small, leading to scale- independent results (plateau). In practice such a condition is hard to satisfy: The 
right inequality is extended to (2 — 5)/a^ leading to lattice artifacts in our results that are of 0{a^p^). Fortunately 
our perturbative calculations allow us to subtract the leading perturbative O(a^) lattice artifacts which alleviates the 
problem. To remove the remaining 0{a?p^) artifacts we extrapolate linearly to a^p^ = as demonstrated in Figs. 
13 -15. The statistical errors are negligible and therefore an estimate of the systematic errors is important. We note 
that, in general, the evaluation of systematic errors is difficult. The largest systematic error comes from the choice 
of the momentum range to use for the extrapolation to a^p^ = 0. One way to estimate this systematic error is to 
vary the momentum range where we perform the fit. Another approach is to fix a range and then eliminate a given 
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momentum in the fit range and refit. The spread of the results about the mean gives an estimate of the systematic 

error. In the final results wc give as systematic error the largest one from using these two procedures which is the 
one obtained by modifying the fit range. We choose the same momentum range in physical units for all /3-values and 
we thus extract all rcnormalization constants using the same physical momentum range, ~ 15 — 32 (GeV)^. This 
momentum range has been chosen so that we are in a region where an approximate plateau is seen at each (3. We 
also note that the 0{a?) perturbative terms which we subtract, decrease as f3 increases, as expected and the values 
extracted from subtracted and unsubtracted data agree when extrapolated to a = 0. The momentum range in lattice 
units at each /3 is as follows: /3 = 3.9 : c?p'^ 3 - 5, /3 = 4.05 : a^p^ - 1.9 - 3, /3 = 4.20 : aV ~ 1-2 - 2.5 and as can 
be seen in Figs. 13 -15 within these ranges the data fall on a straight line of a small slope. 

Our final results for the Z-factors in the MS-schcmc at 2 GeV arc given in Table VII, which have been obtained by 
extrapolating linearly in a?p^, using the fixed momentum range ~ 15 — 32 (GeV)^. 



p 


.^DVl 




Zt> Al 


Zt)A2 


3.90 
4.05 
4.20 


0.970(34) (26) 
1.033(11)(14) 
1.097(4) (6) 


1.061(23)(29) 
1.131(23)(18) 
1.122(7)(10) 


1.126(22)(78) 

1.157(9)(7) 

1.158(7)(7) 


1.076(5)(1) 

1.136(5) 

1.165(5)(10) 



TABLE VII: Rcnormalization constants Zdv and Zua in the MS scheme. The above values have been obtained by extrapolating 
linearly in a^p^ . Statistical errors are are shown in the first parenthesis. The error in the second parenthesis is the systematic 
error due to the extrapolation, namely the difference between results using the fit range p ~ 15 - 32 (GeV)^ and the range 
~ 17 — 24 (GeV)^. An error smaller than the last digit given for the mean value is not quoted. 



VI. CONCLUSIONS 

The values of the renormalization factors for the one-derivative twist-2 operators are calculated non-perturbatively. 
The method of choice is to use a momentum dependent source and extract the renormalization constants for all the 
relevant operators. This leads to a very accurate evaluation of these renormalization factors using a small ensemble 
of gauge configurations. The accuracy of the results allows us to check for any light quark mass dependence. For all 
the renormalization constants studied in this work we do not find any light quark mass dependence within our small 
statistical errors. Therefore it suffices to calculate them at a given quark mass. Wc also show that, despite of using 
lattice spacing smaller than 1 fm, 0{a?) effects are sizable. Wc perform a perturbative subtraction of C(a^) terms. 
This leads to a smoother dependence of the renormalization constants on the momentum values at which they are 
extracted. Residual 0{a^p'^) effects are removed by extrapolating to zero. In this way we can accurately determine 
the renormalization constants in the RI'-MOM scheme. In order to compare with experiment we convert our values 
to the MS scheme at a scale of 2 GeV. The statistical errors arc in general smaller than the systematic. The latter 
are estimated by changing the window of values of the momentum used to extrapolate to a^p^ = 0. Our final values 
are given in Table VII. 

VII. ACKNOWLEDGMENTS 

This work was partly supported by funding received from the Gyprus Research Promotion Foundation under 
contracts EPYAN/0506/08, and TECHNOLOGY/eEniS/0308(BE)/17. 



[1] S. Durr et al.. Science 322, 1224 (2008). 

[2] C. Alexandrou et al. (ETM), Phys. Rev. D80, 114503 (2009), arXiv:0910.2419. 
[3] X.-D. Ji, J. Phys. G24, 1181 (1998), hep-ph/9807358. 

[4] P. Hagler et al. (LHPC), Phys. Rev. D68, 034505 (2003), hepdat/0304018. 
[5] P. Hagler et al. (LHPC), Phys. Rev. D77, 094502 (2008), arXiv:0705.4295. 
[6] S. N. Syritsyn et al. (2009), arXiv:0907.4194. 

[7] D. Brommel et al. (QCDSF-UKQCD), PoS LAT2007, 158 (2007), arXiv:0710.1534. 

[8] C. Alexandrou et al. (ETMC) (2008), arXiv:081 1.0724. 

[9] T. Yamazaki et al., Phys. Rev. D79, 114505 (2009), arXiv:0904.2039. 



25 



[10] C. Alexandrou et al., PoS LAT2009, 145 (2009), arXiv:0910.3309. 
[11] C. Alexandrou et al., PoS LAT2009, 136 (2009). 

[12] M. Constantinou, H. Panagopoulos, and F. Stylianou, PoS LAT2009, 205 (2009). 

[13] C. Alexandrou, M. Constantinou, T. Korzec, H. Panagopoulos, and F. Stylianou, in preparation. 

[14] P. Wcisz, Nucl. Phys. B212, 1 (1983). 

[15] M. Constantinou, V. Lubicz, H. Panagopoulos, and F. Stylianou, JHEP 10, 064 (2009), arXiv:0907.0381. 

[16] G. P. Lepage and P. B. Mackenzie, Phys. Rev. D48, 2250 (1993), hep-lat/9209022. 

[17] P. Dimopoulos et al., PoS LAT2007, 241 (2007), arXiv:0710.0975. 

[18] P. O. L. Zhaofeng, V. Morcnas, private communication. 

[19] M. Gockeler ct al., Nucl. Phys. B544, 699 (1999), hcp-lat/9807044. 

[20] P. de Forcrand, Nucl. Phys. Proc. Suppl. 9, 516 (1989). 

[21] M. Constantinou et al., JHEP 08, 068 (2010), 1004.1115. 

[22] J. A. Gracey, Nucl. Phys. B667, 242 (2003), hep-ph/0306163. 

[23] E. G. Floratos, D. A. Ross, and C. T. Sachrajda, Nucl. Phys. B129, 66 (1977), errarum, ibid., B139, 545-546 (1978). 



